function poisson_fe(n)
%build our grid of nodes that we need to find
if nargin < 1
  n = 10;
end

g=zeros(n,n);

g=sym(g);

%build symbolic matrix by name
for ni = 1:n
    for nj = 1:n
       str=sprintf('z%d%d', ni, nj);
       g(ni,nj)=sym(str);
    end
end


x=linspace(0,4,n);
y=linspace(0,4,n);

g(1,:) = zeros(1,n);
g(:,1) = zeros(n,1);
g(n,:) = x;
g(:,n) = y;


syms S;
for xi= 1:n-1
  for yi=1:n-1
      %upper and lower triangle vertices
      Tl=[1 x(xi) y(yi+1); 1 x(xi) y(yi); 1 x(xi+1) y(yi)];
      Tu=[1 x(xi) y(yi+1); 1 x(xi+1) y(yi+1); 1 x(xi+1) y(yi)];
      
      %coeficient matrices
      coef1l = Tl\[1 0 0]';
      coef1u = Tu\[1 0 0]';
      coef2l = Tl\[0 1 0]';
      coef2u = Tu\[0 1 0]';
      coef3l = Tl\[0 0 1]';
      coef3u = Tu\[0 0 1]';
      
      
      %Calculate lower triangle
      S=S+(g(xi,yi+1)*coef1l(2)+g(xi,yi)*coef2l(2)+g(xi+1,yi)*coef3l(2))^2 + ...
        (g(xi,yi+1)*coef1l(3)+g(xi,yi)*coef2l(3)+g(xi+1,yi)*coef3l(3))^2;
    
      %Calculate upper triangle
      S=S+(g(xi,yi+1)*coef1u(2)+g(xi,yi)*coef2u(2)+g(xi+1,yi)*coef3u(2))^2 + ...
        (g(xi,yi+1)*coef1u(3)+g(xi+1,yi+1)*coef2u(3)+g(xi+1,yi)*coef3u(3))^2;
    
  end
end


%Build the system of equations for minimizations as a cell array
k = 0;
for i=2:n-1
  for j=2:n-1
    k=k+1;
    minimizer{k} = diff(S,g(i,j));
  end
end



%Pass the cell array in so we can solve
sol = solve(minimizer{:});

%get names and values
solnames = fieldnames(sol);
sol = struct2cell(sol);

%substitute in the names and values
zs=g;
for k=1:size(sol)
    zs=subs(zs,solnames{k},sol{k});
end
zs
z=double(zs);

for xi = 1:n-1
  for yi = 1:n-1
     patch([x(xi); x(xi); x(xi+1)],[y(yi); y(yi+1); y(yi)], ...
          [z(xi,yi); z(xi,yi+1); z(xi+1,yi)], ...
          [z(xi,yi); z(xi,yi+1); z(xi+1,yi)]);
     patch([x(xi+1); x(xi); x(xi+1)],[y(yi+1); y(yi+1); y(yi)], ...
          [z(xi+1,yi+1); z(xi,yi+1); z(xi+1,yi)], ...
          [z(xi+1,yi+1); z(xi,yi+1); z(xi+1,yi)]);
      
  end
end
grid on;
view(20,20);
